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Abstract 

The origm of spurious solutions in computational electromagnetics, which violate 
the divergence equations, is deeply rooted in a misconception about the first-order 
Maxwell’s equations and in an incorrect derivation and use of the curl-curl equa- 
tions. The divergence equations must be always included in the first-order MeixweU’s 
equations to maintain the ellipticity of the system in the space domain and to guar- 
antee the uniqueness of the solution and/or the accuracy of the numerical solutions. 
The div-curl method and the least-squares method provide rigorous derivation of 
the equivalent second-order Maxwell’s equations and their boundary conditions. The 
node-based least-squares finite element method(LSFEM) is recommended for solving 
the first-order full Maxwell equations directly. Examples of the numerical solutions 
by LSFEM for time-hajmonic problems are given to demonstrate that the LSFEM is 
free of spurious solutions. 
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1 Introduction 

The occurrence of spurious solutions in computational electromagnetics has been 
known for more than two decades, and elimination of such non-physical solutions is 
still a subject of great interest. The noted feature of these fictitious solutions has been 
their violating the divergence-free conditions in cases where the physical solution is 
completely solenoidal. There is a vast body of reports about spurious solutions asso- 
ciated with the finite element method, see e.g., Cendes and Silvester [10], Bird [3], 
Ikeuchi et al. [22], Davies et al. [14], Rahman and Davies [52] [53], Winkler and Davies 
[68], Webb [66], Welt and Webb [67], Koshiba et al. [30] [31], Ise et al. [21], Rahman 
et al. [54] and Schroeder and Wolff [56]. The majority of spurious solutions has been 
found in the context of eigenvalue analysis. A spurious mode does not correspond 
to the physical modes which the waveguide or resonator imder consideration actu- 
ally supports. The spurious mode problem is severe and often renders the numerical 
solution useless. The spurious solutions have been also revealed in boundary- value 
problems, see, e.g., Crowley et al. [13], Pinchuk et al. [50], Wong and Cendes [69] [70] 
and Paulsen and Lynch [49]. 

The phenomenon of spurious solutions is not exclusive with the finite element 
method. This phenomenon has been edso reported in the context of the fimte differ- 
ence method, see e.g.. Con and Davies [12], Mabaya et el. [36], Schwieg and Bridges 
[57] and Su [61], the boundary element method, see e.g., Ganguly and Spielman [17] 
and Swaminathan et al. [62], and the spectral method, see Fanar and Adams [15]. 
This fact itself undermines the common belief that the spurious solution is a result of 
numerical process. In our opinion, the trouble of spurious solutions in computational 
electromagnetics is deeply rooted in a misconception of the first-order Maxwell’s equa- 
tions and in an incorrect derivation and use of the second-order curl-curl equations. 
We s^ree with Mur [44] [45] that spurious solutions can only be avoided by a correct 
formulation of the problem to be solved. 

In terms of the type of differential equations to be solved, conventional numerical 
methods in computational electromagnetics may be classified into four categories: (1) 
those based on the first-order curl equations; (2) those based on the second-order 
curl-curl equations; (3) those based on the Helmholtz equations; (4) those based on 
the potentials. 

The most widely used numerical method for the solution of time-dependent elec- 
tromagnetic problems has been the finite-difference time-domain (FD-TD) scheme 
developed by Yee [72] and extensively utilized and refined by Taflove and Umashankar 
[63] and Kunz and Luebbers [33], as well as others. In the Yee scheme, only the two 
Maxwell’s curl equations are solved. Some other time-domain methods are also based 
on the two Maxwell’s curl equations, such as the finite volume method developed by 
Shankar et al. [58], the finite difference and finite volume methods by Shang [59] and 
Shang and Gaitonde [60], and the finite element methods by Mei and his colleagues 
[8], Madsen and his colleagues [37] [34], Noack and Anderson [47] and Ambrosiano et 
al. [1]. In general, these approaches do not produce noticeable spurious solutions. This 
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is attributed to the fact that by taking the divergence of the Faraday and Ampere 
laws, one finds that these divergence-free conditions will be satisfied for all time if 
they are satisfied initizdly. However, it is not so easy to satisfy them initially in these 
methods. In fact, in these papers the satisfaction of divergence-free conditions was 
not even considered except by Shang and Gaitonde [60] who seriously examined the value 
of divergence of the computed magnetic field. 

In the original full Maxwell’s equations, when the constitutive relations are speci- 
fied, for three dimensional cases there are eight first-order equations but only six un- 
known vector components, and for two dimensional TE and TM cases four equations 
and three unknowns. That is, the number of equations is larger than the number 
of the unknown functions. For this re<ison, it is traditionally believed that the full 
first-order Maxwell’s equations are “overdetermined” or “overspeafied”, and the two 
divergence equations are thus regarded as “auxihary” or “dependent” and are often 
ne^ected in numerical computation. 

The first-order full Maxwell’s equations have a mathematical structure in which 
the fundamental ingredient is the div-curl system that looks “overdetermined”. A 
similar situation exists in fluid dynamics, see Jiang et d. [27]. By introducing a dummy 
variable( Chang and Gunzburger [11]), however, it can be shown that the div-curl 
system is not “overdetermined”. In this paper, we use this technique to study the full 
M a-rw ^ll’s equations and show that they are properly determined, that is, the two 
divergence equations should not be ignored regardless in dther the static or in the 
time- varying cases. 

In dectromagnetics, there are mainly two reasons why the second-order curl- curl 
equations are often employed. First, it is hard for conventional numerical methods to 
deal with the non-self-adjoint first-order derivatives. Second, in the curl-<nirl equations 
the electric fidd and the magnetic field are decoupled. The curl-curl equations are 
derived from the first-order Maxwell’s curl equations by applying the curl operator. 
It seems that no one has addressed a very important issue: the curl-curl equations 
obtained by simple differentiation without additional equations and boundary condi- 
tions admit more solutions than do its progenitors. In order to derive an equivalent 
higher-order system from a system of vector partial differentid equations, one should 
use the div-curl method that is based on the theorem: if a vector is divergence-free 
and curl-free in a domain, and its normal component or tangential components on 
the boundary is zero, then this vector is identically zero. In other words, the curl and 
the divergence operators must act together with appropriate boundary conditions to 
guarantee that there are no spurious solutions in the resulting higher-order equations. 
In this paper, this div-curl method ori^naUy developed by Jiang et al. [27] is em- 
ployed to derive the second-order system of time-dependent Maxwell’s equations and 
its boundary conditions, and to show that the divergence equations and additional 
boundary conditions must be supplemented to the curl-curl equations. 

The common approach to removing spurious vector modes in the curl-curl equa- 
tions is to modify the variational functional by penalizing the non-zero divergence. 
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Tlie key to success with this so-called penalty method, first used by Hata et al [20] 
and Rahman and Davies [53], depends on the choice of the correct penalty factor — 
values too small or too large do not eliminate spurious solutions. Unfortunately, this is 
an ad hoc and problem- dependent treatment and there has been a lack of systematic 
study of the rationale for selecting this parameter for general problems. 

Recently, the edge element method of Nedelec [46], see e.g., Bossavit and Verite 
[5], Hano [19], Mur and Hoop [43], Barton and Cendes [2], Bossavit [4], Bossavit 
and Mayergoyz [6], Monk [41], Jin [28], Volakis et al. [65] and the references therein, 
li3U! been advocated, because it is believed to be a cure for many difficulties that 
are encountered when attempting to solve electromagnetic field problems by using 
conventional node-based finite elements. Apart &om the fact that such an approach 
can only be used in the simple divergence-free case, edge elements violate the normal 
field continuity between adjacent elements in the homogeneous material domain, see 
Mur [45] for comments and an example. The accuracy of edge elements is lower than 
that of the nodal elements for the same number of unknowns, or the computational 
cost of edge elements is much higher than that of nodal elements for the same accuracy, 
see Mur [45] and Monk [42]. The edge element method also needs non-conventional 
meshing and postprocessing which are not normally available. Moreover, Ross et al. 
[55] reported that the edge element method broke down for large-scale computations 
due to the fact that edge elements c 2 mnot represent purely TE fields. 

It is well known that the solution of the Helmholtz equations with proper boundary 
conditions is free of spurious modes, see Mayergoyz and D’Angelo [38]. The key issue in 
the Helmholtz method is how to specify proper boundary conditions. In this paper, 
we use the div-curl method and the least-squares method to derive the Helmholtz 
equations and their boundary conditions, and show that the divergence equations 
need to be enforced only on a part of boundary, ^lnd they will be implicitly satisfied 
in the domain. We also give a Galerkin variational formulation which corresponds to 
the Helmholtz equations. This theoretically justifies that the penalty parameter s in 
the penalty method should be equal to one. 

The potential approach is widely used in computation of static fields and eddy 
currents. Although the potential approach, see e.g., Boyse et al. [7] for time-harmonic 
problems, does not pve rise to spurious modes, it involves difficulties related to the 
appropriate gau^ng method and the loss of accuracy of the calculated field intensity 
from the potentials by the numerical differentiation. 

This paper emphetsizes that in any case the divergence equations must be included 
explicitly or implicitly as a part of the formulation for electromagnetic problems. 
However, it is not so easy to combine the divergence equations in conventional meth- 
ods. Attempts to satisfy the divergence-free equations by using edge elements merely 
complicate the situation by introducing the need to impose an additional condition 
of normal field continuity. 

This paper shows that the satisfaction of the divergence equations and the elimina- 
tion of spurious solutions can be achieved easily by the application of the node-based 
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least-squaxes finite element method (LSFEM). We believe that the LSFEM is the 
best choice among the available methods for numerical solution of many problems 
in electromagnetics, since it is simple, universal, optimal, robust and eifident. The 
LSFEM is based on the minimization of the residuals in first-order partial differen- 
tial equations. The LSFEM has been successfully applied to various fluid dynamics 
problems, see e.g., Jiang et al. [24] [26], Tang and Tsang [64] and Lefebvre et al. 
[35]. The LSFEM is naturally suitable for the first-order full Maxwell’s equations. 
The preliminary results of LSFEM for time-domain scattering wave problems can be 
found in Wu and Jiang [71]. The theory and the least-squares method for the div-curl 
system discussed in this paper can be employed to directly solve static electric or 
magnetic fields without introdudng the potentials and gauging. In the last section 
of this paper we briefly discuss the general formulation of the LSFEM and apply it 
to time-harmonic problems. Numerical examples are given to demonstrate that the 
LSFEM is free of spurious solutions. 


2 The Div-Curl System 


In this section, we study the div-curl system. We shall show that the three di- 
mensional div-curl system is not “over determined”. We shall introduce the div-curl 
method to derive a second-order system equivalent to the div-curl system. We shall 
show why the least-squares method is the best method for the solution of the div-curl 
system. The technique and the procedure developed here will be applied to dealing 
with the Maxwell’s equations. Since the static Maxwell’s equations are typical div- 
curl systems, the least-squares method introduced in this section can be applied to 
the direct solution of static electric or magnetic fields. 


2.1 Basic Theorems 


First we introduce some notations which axe common in functional analysis. These 
notations will help us to write the mathematical formulations more concisely. Let 
II C K® be a bounded, simply connected, convex and open domain with a piece- 
wise smooth boundary F = Fi U F 2 . Either Fi or F 2 , not both, may be empty. Also 
Fi and F 2 must have at least one commom point, x = (x,y,z) be a point in H, n 
be a unit outward normal vector and r be a tangential vector to F at a boundary 
point, respectively. L 2 { 0 ,) denotes the space of squaxe-integrable functions defined on 
II equipped with the inner product 



uvdSl 


and the norm 


l^llo.n = («>“)• 



6 


H'{Cl) denotes the Sobolev space of functions vrith square-integrable derivatives of 
order np to r. || • ||r,n denotes the usual norm for For vector-valued functions 

u with m components, we have the product spaces 


with the inner product 


and the corresponding norm 


(u,v) = /^» 


• vdft 


|u|lo.n = E ll“illS.n. Il«!lr.n = E ll«illr.n- 
j=i i=i 


Further we define 


<u,v >r= 




uvdT. 


When there is no chance for confusion, we will often omit the measure fl or F from 
the inner product and norm designation. 

Throughout the paper C denotes a positive constant dependent on ft with possibly 
different values in each appearance. 

The following theorems are essential in this paper. 


Theorem 1. If u € fl^^(ft)®, then n x u = 0 on F2^0-<=^ii‘V xu — Oon F 2 . 

Here the notation stands for “leading to and vice versa”. The proof of The- 
orem 1 is straightforward by using the Stokes theorem, see Pironneau [51] or Jiang et 
al. [27]. 


Theorem 2 (Friedrichs’ Div-Curl Inequality). Every function u of H^(ft)® with 
n • u = 0 on Fi and n x u = 0 on F 2 satisfies: 

Hull? <C(1|V -0115 + IIVxullS), (2.1) 

where the constant C > 0 depends only on ft. 

The proof of Theorem 2 involves lots of mathematics. We refer to Girault and 
Raviart [18], Krizek and Neittaanmaki [32] and Jiang et al. [27]. This theorem implies 
that the div-curl norm appearing in the right-hand side of (2.1) is equivalent to the 
norm. This theorem plays a key role in the analysis of the least-squares method. 
From Theorem 2, we can immediately obtain the following theorem which is the basis 
of the div-curl method for deriving higher-order vector equations: 

Theorem 3 (The Div-Curl Theorem). If u € H^(ft)® satisfies 

V • u = 0 in ft, 
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then 


V X u = 0 
n • u = 0 
n X u = 0 

u = 0 


in r2, 
on Fi, 
on F2, 

m Him 


This theorem can also be proved easily by introducing the potential. 

Theorem 4 (The Gradient Theorem). K 5 G H^{U) satisfies 

'Vg = 0 in Q, 

g = 0 on Fi ^ 0 (or on F2 / 0), 

then 


^ = 0 in 12. 


The validation of Theorem 4 is obvious. In fact, ^ = 0 needs to be specified only 
at any point in the domain or on the boundary. This theorem will be used to derive 
the higher-order equations which are equivalent to a scalar equation. 


2.2 The Div-Curl System 

Let us consider the following three-dimensional div-curl system; 


V X u = w in ft. 

(2.2o) 

V • u = p in ft. 

(2.26) 

n • u = 0 on Fi, 

(2.2c) 

n X u = 0 on F2, 

(2.2d) 


where the given vector function u G i/2(fl)* must satisfy the following compatibility 
conditions: 


V • u> = 0 tn 12, 
n • w = 0 on F2, 

J n • ufds = 0. 


(2.3a) 

(2.36) 

(2.3c) 


If F2 is empty, then the given 
condition: 


scalar function 



= 0 . 


L 2 (Cl) must satisfy the compatibility 

(2.3d) 
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At first glance, System (2.2) seems “overdetermined” or “overspedfied”, since 
there are four equations and three unknowns. For this reason, indeed, solving (2.2) 
is not trivial by conventional finite difference or finite element methods. However, 
after careful investigation we shall find that System (2.2) is properly determined and 
elliptic. 

By introducing a dummy variable i?, System (2.2) cjin be written as 


Vi? + V X u = w in Cl, (2.4o) 

V • u = /) in Cl, (2.46) 

n-u = 0 on Fi, (2.4c) 

1 ? = 0 on Fi, (2.4d) 

n X u = 0 on Fj. (2.4e) 


Notice that we impose i? = 0 on Fi, and do not spedfy any boimdary condition for 
the dummy variable i? on F 2 . 

By virtue of Theorem 3, Eq. (2.4a) is equivalent to the following equations and 
boundary conditions: 


V X (Vi? + Vxu — w) = 0 

in fl. 

(2.5a) 

V • (Vi? + Vxu — a>) = 0 

tn fl. 

(2.56) 

n X (Vi? + V XU — u;) = 0 

on Fi, 

(2.5c) 

n • (Vi? + Vxu — co)=0 

on F 2 . 

(2.5d) 


Taking into account the compatibility conditions (2.3a) and (2.3b), the boundary 
condition (2.4e) and Theorem 1, Eq. (2.5b), (2.4d) and (2.5d) lead to 

Ai? = 0 in Cl, (2.6a) 


1 ? = 0 on Fi, 


1 ^- 


on I 2 . 


( 2 . 66 ) 

(2.6c) 


From (2.6) we know that 1 ? = 0 in fl. That is, the introduction of 1 ? into (2.2) does 
not change anything, and thus System (2.4) with four equations and four unknowns 
is indeed equivalent to System (2.2). 

Now let us classify System (2.4). In Cartesian coordinates the equations in System 
(2.4) are given as 

dw dv _ 
dx dy dz *’ 


d‘& du dw 
dy dz 




dx 


yt 


(2.7) 
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d"d dv du _ 

dz dz dy *’ 

du dv dw _ 

dx ^ dy dz ^ 

We may write System (2.7) in the standard matrix form: 


5u 5u . 

Ai-^ + AoU — f, 


in which 


(0 

0 


0 

1 \ 


' 0 

0 

1 


0 

0 


-1 

0 


0 

0 

0 

1 

0 

1 


0 

0 

, A2 = 

-1 

0 

0 

0 

Vi 

0 


0 

0 / 


. 0 

1 

0 

0 / 


- 

-1 

0 

O' 


/o 

0 

0 

0 \ 

1 


0 

0 

0 

, Ao = 

0 

0 

0 

0 

O 1 


0 

0 

1 

0 

0 

0 

0 

\0 


0 

1 

0 > 


\0 

0 

0 

0 / 



f 



f = 

Wy 

, u = 

V 


Wz 


w 


\ p ) 


Uy 


The characteristic polynomial associated with System (2.7) is 


(2.8) 


det(Ai^ + A.2T] + As^) — det 


c 

-V 


-C V C\ 
0 -i V 
^ 0 C 

V C oj 


ie+v^+cy^o 


for all nonzero real triplets (^,J7,C)> System (2.4) is thus elliptic and properly deter- 
mined. 

The fist-order elliptic system (2.4) has four equations and four unknowns, so two 
boundary conditions on each boundary are needed to make System (2.4) well-posed. 
Here ^ = 0 and n • u = 0 serve as two boundary conditions on Fi; while n x u = 0 
implies that two tangential components of u are zero on F 2 . 

Since System (2.2) is equivalent to System (2.4), and System (2.4) is elliptic and 
properly determined, so is System (2.2). 

Itemark In fact, the compatibility conditions (2.3a, b) can be obtained by appl 3 dng 
the div-curl method to the equation (2.2a). 
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2.3 The Div-Curl Method 

Let us derive a higher-order system which is equivalent to the div-curl system (2.2). 
By virtue of Theorem 3, System (2.2) is equivalent to the following system: 


V X (V X u — w) = 0 in fl, (2.9a) 

V • (V X u — a») = 0 in f2, (2.96) 

nx(Vxu — u;) = 0 on Fi, (2.9c) 

n*(Vxu — <o) = 0 on F 2 . (2.9d) 

V • u = p in n, (2.9e) 

n • u = 0 on Fi, (2.9/) 

n X u = 0 on F 2 . (2.9^) 

Due to the compatibility conditions (2.3a, b), the boundary condition (2.9^) and The- 
orem 1, (2.9b) and (2.9d) £ire satisfied. Therefore, System (2.9) can be simplified 
as 

V X (V X u) = V X CO in fl, (2.10a) 

V • u = /) in fl, (2.106) 

n • u = 0 on Fi, (2.10c) 

nx(Vxu) = nxco on Fi, (2.10d) 

n X u = 0 on F 2 . (2.10e) 


Now at least one thing is made clear by the div-curl method. That is, the curl-curl 
equation (2.10a) cannot stand alone; it must go with the divergence equation (2.10b) 
and the additional Neumann boundary condition (2.10d). 

System (2.10) can be further simplified. By virtue of Theorem 4, Eq. (2.10b) is 
equivalent to the following system of equations and boundary condition (assuming 
that F 2 ^ 0): 

V(V*u— p) = 0 in Cl, (2.11a) 

V • u = p on F 2 . (2.116) 

Taking into account (2.11) and the following vector identity: 

V X V X u = V(V • u) - Au, (2.12) 

System (2.10) can be reduced as 


Au = —V X u; Vp in il, 
V(V • tt — p) = 0 in Cl, 


(2.13a) 

(2.136) 
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n • u = 0 on Fi, 

(2.13c) 

n X (V X u) = n X w on Fi, 

(2.13d) 

n X u = 0 on Fa, 

(2.13e) 

V • u = p on Fa. 

(2.13/) 


The solution of the derived second-order system (2.10) or (2.13) is completely identical 
to the solution of the original div-curl system (2.2), therefore no spurious solution will 
be produced by the system (2.10) or (2.13). Moreover, the divergence equation (2.13b) 
in System (2.13) can be deleted. That is, the divergence equation is implicitly satisfied 
by the equation (2.13a) and boundary conditions (2.13c-f). The rigorous proof of this 
statement will be given by using the least-squares method in the next section. Here 
we give a simple explanation adopted from Mayergoyz and D’Angelo [38]. Let us 


consider a slightly different problem: 

Au = —V X w -f Vp in D, (2.14a) 

n • u = 0 on Fi, (2.146) 

n X (V X u) = n X to on Fi, (2.14c) 

n X u = 0 on Fj, (2.14d) 

V-u-p = 0onF. (2.14e) 

That is, we let the divergence equation be satisfied on the whole boundary. Although 
this condition needs to be specified only on Fj, it is not wrong for it to be enforced on 
F. By taMng the divergence of (2.14a) we obtain a Poisson equation of <^ = V • u - p: 

A<^ = 0 in a (2.15) 


Since ^ = 0 on the whole boundary, <f> must be equal to zero in the whole domain, 
i.e., the divergence equation is implicitly satisfied in the system (2.14). 

2.4 The Least-Squares Method 

Let us introduce a more powerful and systematic method, the least-squares method, 
to solve System (2.2) and to derive a higher-order system without spurious solutions. 
We construct the following quadratic frmctional: 

I:H — 

/(u) = ||Vxu-w||§-l-l|V.u-pl|g, 

where Ti = {u e H^(fi)®|n u = 0onFi,nxu = 0on Fa}. We note that the 
introduction of a dummy variable in Section 2.2 is only for the verification of the 
determination, and it is not required in the least-squares functional J. Taking the 
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variation of I witH lespect to u, and letting Su — v and SI — 0, we obtain a least- 
squares variational formulation of the following type: find u € H such that 

B(u,v) = £(v) Vv € “W, (2.16) 

where B{', •) is a bilinear form of the type 

jB(u, v) = ( V X u, V X v) -I- (V • u, V • v), 

and X(') is a linear form of the type 

X(v) = (w, V X v) -H (p, V • v). 

Obviously the bilinear form B(u, v) is symmetric and continuous. The coerdveness 
of B(u, v) is due to Theorem 2. Therefore, we immediately have 

^\\n\\l < B(u,u) = i(u) < l|u|li(||«||o + IHIo). 

By virtue of the Lax-Milgram theorem, see e.g., Oden and Reddy [48] or Johnson 
[29], in fact we have proved the following theorem. 

Theorem 5. The solution of (2.16) uniquely exists and satisfies: 

||u|[i < C-dloillo + IHIo). (2.17) 


The following theorem about the error estimate is also a direct consequence of the 
above results. 

Theorem 6. The LSFEM based on (2.16) has an optimal rate of convergence and 
an optimal satisfaction of the divergence condition: 

||u-Uh||o<C/i‘+^||u||2, (2.18a) 

II V • (ufc - p)||o < (7h*||u||2, (2.186) 

where Uh is the finite element solution, k is the order of complete polynomiab used 
in the finite element interpolation. 

The error estimate (2.18a) is not totally new. The early results were obtained by 
Fix and Rose [16] for the case F 2 = 0 and by Krizek and Neittaanmaki [32] for the 
case Fi = 0. The two-dimensional numerical results obtained by using the primitive 
variables insteady of the potential and a prehminary analysis can also be found in 
Jiang and Cheii [23] and Carey and Jiang [9]. 

The advantages of the LSFEM over the potential method for solving the div-curl 
system is obvious: the trouble of selecting a proper gauging method is avoided; the 
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electric or magnetic fields are obtained directly without numerical differentiation and 
thus have a higher accuracy; eind the electric or magnetic fields are continuous across 
the element boundaries. 

In order to further understand the least-squares method, we derive the Euler- 
Lagrange equations associated with the least-squares variational formulation (2.16) 
which can be rewritten as: find u G 'H, such that 

(V X u — w, V X v) (V • u — /), V • v) = 0 Vv G “H. (2.19) 

Suppose that u, eind p are sufficiently smooth. By using Green’s formulae, Eq. 


(2.19) can be written as 

(V X (V X u - w), v)-|- < (V X u — a>), n x v >r + 

(— V(V • u — p), v)-f- < (V • u — p),n • V >r= 0 Vv G "H- (2.20) 

Taking into account (2.12) and that v satisfies n • v = 0 on Ti and n x v = 0 on Tj, 

from (2.20) we obtain 

(-Au - V X w -F Vp,v) 

— <nx(V XU — w),v >r, -f- < (V • u — p), n • v >F2= 0 (2.21) 

for all admissible v G W, hence we have the Euler-Lagrange equation and boundary 

conditions: 

Au = —V X u> -h Vp in Q, (2.22a) 

n • u = 0 on Fi, (2.226) 

nx(Vxu) = nxw on Fi, (2.22c) 

n X u = 0 on F 2 , (2.22d) 

V • u = p on F 2 . (2.22e) 


We note that in System (2.22) the divergence equation does not appear in the 
domain. In fact, we have rigorously proved that the solution of the Helmholtz-type 
system (2.22a) under additional boundary conditions (2.22c, e) automatically satisfies 
the divergence equation. We also remark that if F 2 is empty, the divergence equation 
does not even appear on the boundary. The attraction of using the higher-order 
system (2.22) now becomes apparent: one avoids dealing with the divergence condition 
(2.2b) which is implicitly satisfied; instead, one deals with three Poisson equations 
that everyone would rather solve. However, we should mention that if one chooses 
to use the finite difference method to solve (2.22a), the additional natural boundary 
conditions (2.22c,e) must be supplemented. 

Now it is clear that the following four formulations are equivalent to each other: 
(1) the first-order div-curl system (2.2); (2) the least-squares variational formula- 
tion (2.16); (3) The Helmholtz-type system (2.22); and (4) the Galerkin formulation 
(2.21). It turns out that the least-squares method (2.16) for the div-curl equations 
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(2.2) corresponds to using the Galerkin method (2.21) to solve System (2.22) which 
consists of three independent second-order Poisson equations (2.22a) and three cou- 
pled boundary conditions on each boundary, where the original first-order equations 
(2.22c) and (2.22e) serve as the natural boundary conditions, and (2.12b) and (2.22d) 
as the essential boundary conditions. 

Obviously, the least-squares problem is formally equivalent to a higher-order 
problem with additional natural boundary conditions provided by the original first- 
order differential equations. The least-squares method (2.16) is the simplest approach 
among these equivalent methods, because it does not need any additional bound- 
ary conditions. The trial function u and the test function v need to satisfy only the 
oripnal essential boundary conditions. This is one of the reasons why we strongly 
recommend the least-squares method. 

Now we have shown that the three-dimensional div-curl system can have three 
equivalent differential forms: (1) the first-order system (2.2); (2) the curl-curl equa- 
tion (2.10a) which must be accompanied by the divergence equation (2.10b) and 
the additional Neumman boundary condition (2.10d); (3) three uncoupled Poisson 
equations (2.22a) with additional Neumman boundary conditions (2.10c) and (2.10e) 
provided by the original first-order system. 

In the following sections, we will show that Maxwell’s equations have similar 
equivalent forms. 


3 The First-Order Maxwell’s Equations 

In this section we shall show that the first-order full Maxwell’s Equations are not 
“overdetermined”, emd thus the divergence equations should not be ignored. 

3.1 The Basic Equations 

For general time- varying fields, the original first-order full Maxwell’s equations can 
be written as 


V X E -f- = — K*’"*’ in II, (Faraday's Law) (3.1a) 

dt 

V X H — - <tE = in II, (Maxwell — Ampere's Law) (3.16) 

dt 

V • (eE) = />*”’** in 12, (Gauss's Law — Electric) (3.1c) 

V • (pH) = 0 tn 12, (Gauss's Law — Magnetic) (3.1d) 


where E and H are the electric and magnetic field intensities respectively, is the 
imposed source of electric charge density, and J*™** and K*”*** are imposed sources of 
electric and magnetic current density. All imposed sources are pven functions of the 
space and time coordinates. 
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In System (3.1) we have already made use of the following constitutive relations: 

D = cE, 

J = <tE, 

where D is the electric flux density, B is the magnetic flux density and J is the 
electric (eddy) current density; and the constitutive parameters e, fi and a denote, 
respectively, the permittivity, permeability and conductivity of the medium. These 
parameters are tensors for 2 inisotropic media, and may be functions of position and 
time, and may depend on the field intensities. For simplicity, we consider isotropic 
and homogeneous media, therefore they are constant scalars. 

The field equations are supplemented by the boundary conditions: 


n X E = 0 

on Fi, 

(3.1e) 

n . (^H) = 0 

on Fi, 

(3.1/) 

n X H = 0 

on F 2 , 

(3.1s) 

n • (eE) = 0 

on F 2 , 

(3.U) 


where Fj is an electric wall, and F 2 is a magnetic symmetry wall. Here we consider 
only homogeneous boundary conditions, since inhomogeneous boundary terms can 
always be converted into source terms. 

For transient problems, the initial conditions on E and H should also be provided. 
To allow System (3.1) to have a solution, the source terms must satisfy the fol- 


lowing compatibility conditions: 

V • K*"**’ = 0 in fl, (3.2a) 

n.K‘"*' = 0 onFi, (3.25) 

^n.K‘"“’dT = 0, (3.2c) 

V • J"”*’ -I- = 0 in ft, (3.2d) 

at 

n • = 0 an Fj. (3.2e) 


We remark that the compatibility conditions (3.2a,b,d,e) can be obtained by ap- 
plying the div-curl method to the Maxwell’s curl equations (3.1a, b). 
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3.2 The Determination 



Consider the following system augmented by the variables (p and x* 


ot 

12, 

(3.3a) 

Vx + VxH = in 12, 

(3.36) 

V • (eE) = in 12, 


(3.3c) 

V • (/iH) = 0 in 12, 


(3.3d) 

n X E = 0 on Fi, 


(3.3e) 

X = 0, n-(/*H) = 0 on Fi, 


(3.3/) 

n X H = 0 on F 2 , 


(3.3^) 

^ = 0, n • (cE) = 0 on F 2 , 


(3.3/i) 

We shall prove that p and x (3-3) are dummy wiables, i.e., System (3.3) is 

equivalent to System (3.1). In fact, by virtue of the div-curl theorem, Eq. (3.3a) and 
(3.3b) are equivalent to the following equations: 

V X {V^-F V X E + = 0 

in 12, 

(3.4a) 

V • {Vv? -1- V X E -F -t- K*”***} = 0 

in 12, 

(3.46) 

n . { V X E -h -1- K*'"*’} = 0 

on Fj, 

(3.4c) 

n X {V^ -h V X E -h -1- K*'"*'} = 0 

on F 2 , 

(3.4d) 

Vx{Vx + VxH <rE = 0 

in 12, 

(3.5a) 

V . {Vx + V X H - - o-E -T’"'’} = 0 

in 12, 

(3.56) 

n X {Vx + V X H - - (tE - J'”***} = 0 

on Fi, 

(3.5c) 

n . {Vx + V X H - - o-E - T™**} = 0 

on F 2 , 

(3.5d) 


Taking into eiccount the divergence-free condition (3.3d) and the compatibility con- 
dition (3.2a), from (3.4b) we find that 


Ay? = 0, 


in Q. 


(3.6a) 
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Taking into account the boundary condition (3.3e) and Theorem 1, the boimdary 
condition (3.3f), and the compatibility condition (3.2b), from (3.4c) we obtain 

n • Vy> = 0 on Fi. (3.66) 

From (Z.Zh) we have 

<p = 0 on F 2 . (3.6c) 

System (3.6) implies that y? = 0 in Similarly, we can show that = 0 in 
Therefore, tp and x (3-3) are really dummy variables, and thus System (3.3) is 
equivalent to System (3.1). 

The first-order system (3.3) has eight equations, eight unknowns, and four bound- 
ary conditions on each part of the boundary, and thus is properly determined. It is 
valid for static, transient, and time-harmonic cases. 

In static cases, the time-derivative terms in (3.3a) and (3.3b) disappear, and oE 
is included into the given current density. The system (3.3) becomes two independent 
div-curl systems for the electric field and the magnetic field respectively. In Section 
2.2, we have shown that each div-curl system is elliptic. 

In time-harmonic cases, when the time factor is used and suppressed, the 
time-derivative terms become the zero-order terms, and the system (3.3) becomes 
two coupled div-curl systems. The coupling is through the zero-order terms. The 
principle part, i.e., the first-order derivative terms which classify the system, stiU 
have the div-curl structure, and thus the whole system is elliptic. 

In transient cases, the whole system (3.3) is hyperbolic. However, in time-domain 
numerical methods, the time-derivative terms are discretized by explicit or implicit 
finite differences, hence the time-derivative terms become the zero-order terms in the 
space domain. For each time step, the time-discretized system is still elliptic. 

In summary, in all cases. System (3.3) is properly determined, and is eUiptic in the 
space domain. Since System (3.1) is equivalent to System (3.3), it is indeed properly 
determined and also elliptic in the space domain. 

3.3 The Importance of Divergence Equations 

It is commonly believed that the divergence equations (3.1c) and (3.1d) axe “redun- 
dant” for transient and time-harmonic problems, and thus axe ne^ected in compu- 
tation. This misconception is the true ori^n of spurious or inaccurate solutions in 
computational electromagnetics due to the following reasons: 

(1) The original first-order full Maxwell’s equations reflect the general laws of 
physics. They govern all electromagnetic phenomena, no matter whether the problem 
is static, time-harmonic or transient. But the first-order curl equations (3.1a) and 
(3.1b) cannot govern static cases. For static problems, the divergence equations (3.1c) 
and (3.1b) must be explicitly included in the first-order Maxwell’s equations (3.1). 

(2) By taking the divergence of (3.1a), one can conclude only that 5(V-(/iH))/ dt = 
0, that is, V • ifiH) = F{x). In (3.1a) there is no information about this function. 
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Some literature asserts that, if the divergence of /*H is zero at the beginning, it will be 
identically zero forever. The problem then is how one can set V • (AtH) = 0 initially. 
Let us examine the common practice: letting the imtial field intensities be zero in 
the domain, and the boimdary conditions be correctly given on the boundary. In this 
case, the divergence-free condition is significantly violated near the boundary at the 
first time step of the computation. If someone really can make the divergence be zero 
at the beginning, it is in fact equivalent to adding a divergence-free equation into the 
system. The discussion for the electric field runs along the same line. 

(3) The neglect of the divergence equations destroys the ellipticity of Maxwell’s 
equations in the space domain. In each curl system there are only three(odd) equations 
and three(odd) unknowns that cannot be elliptic in the ordinary sense. In general, the 
umnerical methods based on a non-elliptic system without special treatment cannot 
be proved to possess an optimal rate of convergence. A related investigation can be 
found in Jiang and Povinelli [25]. 

(4) The time marching method is often an effective approach to solve steady-state 
non-linear problems where the material properties depend on the electromagnetic 
fields. The curl equations alone are not adequate for this approach. Neither are the 
curl equations appropriate for solving the scattering of waves excited by a pulse wave. 


4 The Second- Order Maxwell’s Equations 

In this section we shall use the div-curl method to derive the second-order 
Maxwell’s equations eind their boundary conditions, and show that the curl-curl equa- 
tions cannot stand alone. We shall give the Galerkin method corresponding to the 
correct second-order Maxwell’s equations. We shall see that this Galerkin method is 
of the same form as the popular Galerkin/penalty method with the penalty parameter 
s = 1. We shall also give a simple least-squares look-alike method to obtain a correct 
variational formulation which rigorously justifies that s = 1 in the penalty method. 

4.1 The Div-Curl Method 

By virtue of the div-curl theorem. System (3.1) is equivalent to 


V X {V X E -f -f = 0 

at 

in fl, 

(4.1o) 

V . {V X E -H + K*”**’} = 0 

at 

tn fl. 

(4.16) 

n . {V X E -h ^ = 0 

at 

<m Fi, 

(4.1c) 

n X {V X E-H = 0 

on Fj, 

(4.1i) 
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V X {V X H ^ <tE - J*’"*'} = 0 

ot 

in ft. 

(4.1e) 

V-{VxH o-E J‘"*'} = 0 

at 

in ft. 

(4-1/) 

n X {V X H o-E = 0 

at 

on Fi, 

(4.1s) 

n ■ {V X H ^ o-E J*"***} = 0 

<7c 

on F 2 , 

(4.U) 

V • (eE) = in fi, 


(4.1i) 

V • (/tH) = 0 iTi fi, 


(4.1j) 

n X E = 0 on Fi, 


(4.U) 

n • (/iH) = 0 on Fi, 


(4.11) 

n X H = 0 on F 2 , 


(4.1m) 

n • (eE) = 0 on F 2 , 


(4.1n) 

Due to the compatibility conditions (3.2), the divergence conditions (4.1t,y) and the 
boundary conditions (4.1fc — n), we may eliminate Equations (4.16), (4.1c), (4.1/) and 
(4.1h) and rewrite System (4.1) as 

Vx{VxE+^^^^}= VxK*’”*’ 
at 

in ft. 

(4.2o) 

V X {V X H ^ <tE} = V X J"”** 

ot 

in ft, 

(4.26) 

V • (eE) = p*"**’ in ft, 


(4.2c) 

V • (/iH) = 0 in ft, 


(4.2i) 

n X E = 0 on Fi, 


(4.2e) 

n • (/xH) = 0 on Fi, 


(4.2/) 

n X (V X H) = n X J*"*** on F 

1 

(4.2s) 

n X H = 0 on F 2 , 


(4.26) 

n • (eE) = 0 on F 2 , 


(4.2i) 

n X (V X E) = — n X K’"*** on 

T 2 . 

(4.2/) 

System (4.2) is completely equivalent to System (3.1), the validation of (4.2) guar- 
antees the validation of (3.1). Therefore, we can use the curl equations in (3.1) to 
decouple E and H in (4.2) as usual, then we obtain 

V X (V X E) + + <^E) = -V X M 

Qjimp 

dt 

(4.3a) 
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V • (eE) = p*"*** in Q, 
n X E = 0 on Fi, 
n • (eE) = 0 on F 2 , 
n X (V X E) = -n X K’’"*’ on Fj, 
id 

V X (V X H) + + <t)?^ = -(s^ + <r)K<"' + V x J‘”' 

V • (pH) = 0 inQ, 
n • (pH) = 0 on Fi, 
n X (V X H) = n X J*™** on Fi, 
n X H = 0 on F 2 . 


(4.36) 

(4.3c) 

(4.3d) 

(4.3e) 

in (4.4a) 

(4.46) 

(4.4c) 

(4.4d) 

(4.4e) 


We note tliat the curl-curl equations in (4.3) and (4.4) cannot stand alone; they must 
be supplemented by the divergence equations and the additional natural boundary 
conditions. In other words, the curl-curl equations admit more solutions than the 
first-order full system. This is the real reason that the numerical methods based on 
the curl-curl equations will give rise to spurious solutions. 

It is difficult to solve a second-order curl-curl equation (4.3a) with the explicit 
constraint of the first-order divergence equation (4.3b), since the problem has an 
overspecified number of partial differential equations and the first-order equation 
(4.3b) is hard to deal with numerically. We shall look for a simple way. By using 
Theorem 4 and the vector identity (2.12), System (4.3) and System (4.4) can be 
reduced to 


-AE + /.|(^ + <tE) = -VxIC^-M^-(;)Vp‘”' inSi, (4.5a) 


V(V • (eE) - p*”***) = 0 in a, (4.56) 

n X E = 0 on Fi, (4.5c) 

V • (eE) = p‘"‘»’ onFi, (4.5d) 

n • (eE) = 0 on Fa, (4-5e) 

n X (V X E) = -n X K*"*** on Fa, (4.5/) 

and 

-AH + (a| + a)^ = -(e|+a)K‘^+VxJ<^ in 0, (4.6a) 

V(V • (pH)) = 0 in fi, (4.66) 
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n • (^H) = 0 on Fi, 

(4.6c) 

n X (V X H) = n X J*’"*’ on Fi, 

(4.6d) 

n X H = 0 on F 2 , 

(4.6e) 

V • (/xH) = 0 on F 2 . 

(4.6/) 


Due to the reasons pointed out in Section 2.3, EJq. (4.5b) and Eq. (4.6b) can be elim- 
inated. That is, the divergence equations (4.5b) and (4.6b) are redundant, they are 
implicitly satisfied by the Helmholtz-type equations (4.5a) and (4.6a) and the bound- 
ary conditions. Therefore, System (4.5) and System (4.6) can be further simplified 
as 


-AE + = -V X 

tn D, 

(4.7a) 

n X E = 0 on Fi, 


(4.75) 

V • (eE) = on Fi, 


(4.7c) 

n • (eE) = 0 on F 2 , 


(4.7d) 

n X (V X E) = — n X K**"** on F 2 , 


(4.7e) 

and 



-AH + = -(4 + .)K- + V X 

in Q, 

(4.80) 

n • (pH) = 0 on Fi, 


(4.86) 

n X (V X H) = n X J**"** on Fi, 


(4.8c) 

n X H = 0 on F 2 , 


(4.8d) 

V • (pH) = 0 on F 2 . 


(4.8e) 

We note that the divergence conditions are required to be satisfied only on 

a put 


of boundary. We will rigorously prove this in Section 4.3 by using the least-squares 
method. As in Section 2.3 for the div-curl system, one may enforce the divergence 
conditions on the whole boimdary F in (4.7) and (4.8) and show that the divergence 
conditions are satisfied in the domain D. 

The Helmholtz-type equations (4.7a) and (4.8a) can be found in all text books 
on electromagnetics. However, it seems that all these books (except for Mayyergoyz 
and D’Angelo [38] who got the same conclusion as ours for a special case) claim that 
the Helmholtz equation must be solved vnth the divergence equation satisfied in the 
whole domain and do not mention that it needs additional botmdary conditions. Our 
rigorous derivation using the div-curl method shows that the Helmholtz equation can 
stand alone, zuid the divergence equation should be satisfied only on a part of the 
boundary. 

The advantages of using the Helmholtz equation over the curl-curl equation are 
obvious: one avoids the difficulty involving the explicit satisfaction of the divergence 
equations, instead one solves three decoupled second-order equations with coupled 
boundary conditions. 
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4.2 The Galerkin Method 

One may elect to use, for example the finite difference method, to solve the Helmholtz- 
type systems (4.7) or (4.8). Usually, the finite difference method is based on rectangu- 
lar structural grids. In this case, for example, the divergence condition (4.8e) can be 
simplified as the Neumann boundary condition (see Mayergoyz and D’Angelo [38]): 

= 0 on Ta. 
an 

For complex geometry it is not straightforward to implement the Neumann boundary 
condition in the finite difference method. By using the finite element method based 
on a variational principle, even the divergence conditions on the boundary do not 
appear. In the following we derive the variational formulation corresponding to (4.7). 

By taking into account the vector identity (2.12), the Galerkin formulation asso- 
ciated with (4.7) is: find E satisfying (4.7b) and (4.7d) such that 

(V X {V X E -f K™*’}, E*)-l- < V X E -1- x E* >r. 


+(_V{V • E - E*)-|- < V • E - ■ E* >r. 




(4.9) 


for all E* satisfying (4.7b) and (4.7d). By virtue of Green’s formula, the statement 
(4.9) can be simplified to a more symmetric form: find E satisfying (4.7b) and (4.7d) 
such that 


(V X E, V X E’) + (V ■ E, V . E*) -I- + «rE},E*) = 

-(K’’"^ V X E*) + ir^/e, V • E*) - (a*^ E*) (4.10) 

for all E* satisfying (4.7b) and (4.7d). 

For time-harmonic(eigenvalue) problems with tr = 0, the variational formulation 
takes the form 


(V X E, V X E*) -1- (V • E, V • E*) - o;Ve(E,E*) = 0, (4.11) 

where w is the angular frequency. The formulations for the magnetic field are similar. 

The variational formulations (4.10) and (4.11) are of the same structure as the 
most popular Galerkm/penalty formulations in the literature. However, in contrast 
to the commonly used penalty formulation, there is no free parameter in the Galerkin 
formulation (4.10) and (4.11). In other words, the penalty parameter s = 1 should be 
chosen in the penalty method in order for the penalty method to correspond to the 
Helmholtz-type equations (4.7). 
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4.3 The Least-Squares Look-Alike Method 

In Section 4.1 the div-curl method is employed to derive the second-order (Helmholtz- 
type) Maxwell’s equations and their boundary conditions that guarantee no spurious 
solutions. But there we cannot make sure that the divergence conditions should be 
spedlied only on a part of boundary. In this section we give a more powerful method 
to derive equivalent higher-order equations and rigorously prove the statement made 
in Section 4,1. 

Consider the following div-curl system for the electric field: 


5(pH) 

dt 

(4.12a) 

V • E = p’"‘»’/£ 

in fl, 

(4.126) 

n X E = 0 

on Fi, 

(4.12c) 

n • (eE) = 0 

on F 2 , 

(4.12d) 


where H is assumed to be known and to satisfy Eq. (3.1b) and the boundary conditions 
(3.1/) and (3.1p), and the source terms satisfy the compatibility conditions (3.2a-e). 
In other words, when the magnetic field and the sources are given, the solution of 
(4.12) will give the corresponding electric field. Obviously, System (4.12) is a typical 
^v-cnrl system that has been investigated in Section 2. 

Following the steps in Section 2.4, we can derive the variational formulation which 
corresponds to System (4.7). We define the quadratic functional: 

7(E) = II V X E -H + ||V - E - ||", 

in which E satisfies the boundary conditions (4.12c, d). The minimization of I leads 
to the variational formulation: 

( V X E -f -I- V X E*) -f- ( V • E - V • E*) = 0, (4.13) 

at 

where E* = 5E and satisfies the same boundary conditions as E. Since H satisfies 
(3.1b) and (3. Ip), from (4.13) we have 

(V X E, V X E*) -K V . E, V • E*) -I- E*) = 

-(K*’”P,V X E*) ^-(p’"‘V£,V•E•)-(/l^J‘"‘^E*), ' (4.14) 

which is exactly the same as (4.10). By using Green’s formula, from (4.14) we can 
obtain the Euler-Lagrange equation (4.7a) and the natural boundary condition (4.7c) 
and (4.7e). That is, the correctness of (4.7) or (4.8) is completely proved. 
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Now we understand that the variational formulation (4.14), the Helmholtz-type 
equation (4.7a) and its boundary conditions, and the first-order system (4.12) are 
equivalent to each other. However, the finite element method based on (4.14) has 
superior advantages: the divergence condition (4.12b) is automaticaUy satisfied, the 
test and trial functions are reqtdred to satisfy only the essential boundary conditions 
(4.12c,d). 

We remark that the procedure to obtain the formulation (4.14) is not a true least- 
squares approach, because (1) we have assumed that H is given and satisfies (3.1b), 
and hence H is not subject to the variation; (2) the true least-squares method always 
leads to a symmetric bilinear form; here the <r related term is not symmetric. Hven 
so, this procedure is mathematically justifiable. It is nothing but a rigorous method 
to derive the Galerkin variational formulation corresponding to the Helmholtz-type 
equations (4.7a) and their boimdary conditions. All derivation provided in this section 
has rigorously proved that the penalty parameter in the Galerkin/penalty method 
should be equed to one. 

5 The Least-Squares Method for First-Order Maxwell’s 

Equations 

In Section 2.4 we have introduced the least-squares method for the pure div-curl 
system governing static field problems, and in Section 4.3 the least-squares look alike 
method for the div-curl system describing time-dependent single(electric or magnetic) 
field problems, and demonstrated the power of the least-squares method. In this 
section we briefly give the formtdations of the LSFEM for the general first-order 
partial differential equations, and apply it to the solution of the time-harmonic first- 
order Maxwell’s equations. 

5.1 The General Formulation 

The least-squares method for the linear operator equation An = f formally is equiva- 
lent to the solution of the higher-order equation A*Au = A*f with An = f serving as 
an additional natural boundary condition, where A* is the adjoint of A in the inner 
product generated by the norm. When directly applied to second-order equations 
this approach requires the use of finite elements and leads to ill-conditioned dis- 
crete systems. In order to use simple elements and obtain a better- conditioned 
algebraic system, the least-squares method discussed here is based on the first-order 
system. The formulation of the least-squares finite element method for general first- 
order steady-state boundjiry- value problems can be foimd in Jiang and Povinelli [24]. 
This formulation can be directly applied to the solution of the first-order steady-state 
and time-harmonic Maxwell’s equations. For time-dependent problems one always can 
use an appropriate finite difference method in the temporal domain, such as the back- 
ward Euler scheme or the Crank- Nicolson scheme, to discretize the time- derivative 
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terms so that in each time-step the problems are converted into botmdary- value prob- 
lems. For completeness, we briefly derive the general least-squares formulation. 

We consider the linear boundary- value problem: 

Au = f in (5.1a) 

Bu = g on r, (5.16) 

where A is a first-order partial differential operator: 

"■* du 

Au = 53 -A-i Q H A-oU, (5.2) 

in which Q G is a boimded domain with a piecewise smooth boundary F, = 
2 or 3 represents the number of space dimensions, u^ = (ui,« 2 j — '“m) is a vector of m 
unknown functions of x = (®i, ..., x^^), Aj and Aq are n x m matrices which depend 
on X, f is a given vector- valued function, B is a boundary algebraic operator, and 
g is a pven vector-valued function on the boundary. Without loss of generality we 
assume that the vector g is null. We should mention that the number of equations n 
in the system (5.1a) must be greater than or equal to the number of imknowns m. 

Considering the boundary condition of the boundary- value problem, we also define 
the function space 

V = {v G Bv = 0 an T}. (5.3) 

Let us suppose that f G L 2 (fi) and A : V -+ L 2 (Q). For an arbitrary trial function 
V G V, we define the residual function: 

R = Av — f in Q. (5-4) 

In general the residual R is not equal to zero, except v is equal to the exact solution 
u. The squared distance between Av and f will be nonnegative: 

\\K\\l= /(Av-f)^dB>0. (5.5) 

Jfi 

A solution u to the problem (5.1) Cctn thus be interpreted as a member of V that 
minimizes the squared distance between Av and f: 

0 = I|R(U)|I5 < ||R(v)||5 Vv€V. 

The least-squares method consists of seeking a mininoizer of the squared distance 
||Av-f||§ in V. We write the quadratic functional in (5.5) as 

7(v) = II Av - f\\l = (Av - f, Av - f). (5.6) 

A necessMy condition that u G V be a mininaizer of the functional I in (5.6) is that 
its first variation vanish at u for all admissible v. That is. 


liTn ^ J(u -t- rv) — (Av)^(Au — f)dfl = 0 Vv G V . 
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Tims, the least-squares method leads us to the variational boundary-value problem: 
Find u € V such that 

B(u,v) = F(v) Vv E V, (5.7a) 

where 

B(u,v) = (Au, Av), (5.76) 

F(v) = (f,Av). (5.7c) 

In the finite element analysis, we first subdivide the domain as a union of finite 
elements and then introduce an appropriate finite element basis. Let N„ denote the 
number of nodes for one element and denote the element shape functions. If equal- 

order interpolations are employed, that is, for all unknown variables the same finite 

element is used, we can write the expansion in each element 


<(3C) 


Nn 


Z)V’,-(x) 




\ ttm ' j 


(5.8) 


where (ui,« 2 , ...,ttm)i are the nodal values at the jth node, and h denotes the mesh 
parameter. 

Introducing the finite element approximation defined in (5.8) into the variational 
statement (5.7), we have the linear algebraic equations 

KU = F, (5.9) 


where U is the global vector of nodal values. The globzil matrix K is zissembled from 
the element matrices 

Ke = / (AV»i,AV’2,— ,-A-V’wJ^(A^i,A^ 2,— ,AV>jv«)<K^, (5.10) 

Jn, 

in which Q* C is the domain of the eth element, and T denotes the transpose, and 
the vector F is assembled from the element vectors 

Fe= f (AV»i, A^2,-..,AV»jv,)^f<Kl, (5.11) 

Jn, 


in (5.10) and (5.11) 

’*•' 8tb- 

AV'j = ^ V'iAo. (5-12) 

i=i 

Prom the above derivation we can immediately find out or further prove that: 

(1) the least-squares method is universal for all tjrpes of partial differential equa- 
tions, no matter whether they are elliptic, hyperbolic, parabolic or mixed; the only 
requirement is that they have a unique solution, see Mikhlin [40] and Marchuk [39]; 
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(2) the LSFEM leads to a symmetric positive definite algebraic system which 
can be solved efficiently by matrix-free iterative methods, such as the element-by- 
element preconditioned conjugate gradient method, and thus the parallelization and 
large-scale 3D computation is made easy; 

(3) the LSFEM formulation and its coding are general, therefore for a new problem 
one needs only to supply the coefficients, the load vector and the boundary conditions; 

(4) the LSFEM is robust, no special treatments, such as upwinding, staggered 
grids, and operator splitting etc. are needed; the LSFEM leads to a minimization 
problem rather than a saddle-point problem, thus simple equal-order interpolations 
can be employed; 

(5) the LSFEM cein often be proved to have optimal numerical properties including 
an optimal rate of convergence; 

(6) the LSFEM satisfies the divergence conditions in electromagnetics. 


5.2 Time-Harmonic Fields 

For three-dimensional time-harmonic fields, the first-order full Maxwell’s equations 
can be written as 


V X E -1- jufiH = - 

-K*"’*’ tn ft. 

(5.13a) 

V X H — y weE = 

J**"** in ft. 

(5.136) 

o 

II 

> 

in ft. 

(5.13c) 

<1 

II 

o 

in ft. 

(5.13d) 

where the time factor is used and suppressed, u is the given angular frequency 

and not equal to the resonant frequendes of this problem, E and H axe the com- 
plex electric and magnetic field intensities respectively, J*”*** and K*"*** are imposed 
harmonic sources of electric and magnetic current density respectively. All imposed 
sources are given functions of the space coordinates. For simplidty, we consider ho- 
mogeneous isotropic media, i.e., e cuid ft axe constant scalars. The field equations are 
supplemented by the homogeneous boundary conditions: 

n X E = 0 

on Fi, 

(5.13e) 

n-H = 0 

on Fi, 

(5.13/) 

n X H = 0 

on F 2 , 

(5.13<7) 

n • E = 0 

on F 2 . 

(5.136) 


where Fi is an electric wall, and F 2 is a magnetic symmetry wall. 

To allow System (5.13) to have a solution, the source terms cannot be arbitrary, 
they must satisfy the following compatibility conditions: 


V . K”"" = 0 in ft. 


(5.14a) 
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n . K’"*’ = 0 

on Fi, 

(5.145) 


= 0, 

(5.14c) 

V • = 0 

in f2. 

(5.14d) 

n . = 0 

on F 2 , 

(5.14e) 

f n • J’^dT 

= 0. 

(5.14/) 


Jr 

As in the time-varying cases, the compatibility conditions (5.14a, b) and (5.14d,e) 
raTi be derived by applying the div-cnrl method to the curl equations (5.13a, b). The 
compatibility conditions (5.14c) and (5.14f) can be obtained by applying the Gauss 
divergence theorem to (5.14a) and (5.14d), respectively. 

Separating the real and imaginary parts in (5.13) leads to 


V X Er — w/iHi = - 


in fl. 

(5.15a) 

V X Ei -t- (jJflKr = - 


in Q, 

(5.156) 

V X H^ -t- wsEi = 

Jimp 

in fl. 

(5.15c) 

V X Hi — weE, = 

Jimp 

in f2. 

(5.15d) 

V • E, = 0 

in Q, 


(5.15e) 

0 

II 

> 

in fl. 


(5.15/) 

V • H, = 0 

in ilf 


(5.15^) 

V • Hi = 0 

in 0. 


(5.15A) 


Obviously, System (5.15) is elliptic, since its principle part consists of four div-curl 
systems. For the solution of (5.15) the least-squares variational formulation is: find 
u = (E„Ei,H,,Hi) € n such that 

B(u, v) = I(v) Vv = (e:,e;,h;,h:) g w, ( 5 . 16 ) 

where = {u G x x x ff^(ft)®|n xE = Oonri,nH = 

0onri,nxH = 0on T 2 ,n • E = 0 on 12 }, and 5(-,-) is the bilinear form 

5(u, v) = (V X E,. - wftHi, V X E; — ufiH*) 

-f(V X Ei -f w/iHr, V X E* ■+ 

+{V xHr+ weEi, V X H; -I- weE^) 

-h(V X Hi - u;eE„ V x H* - weE*) 
+(V.E„V-E:)-f(V.Ei,V.E*) 
+(V.B„V.E:)-KV.Ei,V.E;), 


(5.17a) 
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and L{-) is the linear form 

x(v) = (-if’”**', V X e; - w^h;) 

+(-if:"‘^VxE: + a;/iH;) 

+(j;:”‘P,VxH; + a;eET) 

+( V X H; - ueE;). (5.176) 

Obviously, the bilinear form B(u, v) in (5.17a) is symmetric and continuous, and 
the linear form L{v) in (5.17b) is continuous. One may prove that if the frequency 
of the exciting source is not equal to the resonant frequencies of this electromagnetic 
system, then the bilinear form ^(u,u) is coercive. By virtue of the Lax-MUgram the- 
orem, the least-squares solution uniquely exists and the corresponding finite element 
solution is of an optimal rate of convergence. In fact, the following statement is the 
consequence of the coerdveness of the bilinear form B(u, u). We will prove it in our 
future reports. 

The LSFEM based on (5.16) has an optimal rate of convergence and an optimal 
satisfaction of divergence-free conditions: 


||E, - E^fclio < 

||V.E,4o<C'/t‘, 

(5.18a) 

||Ei-Eifc||o<C7h‘+i, 

\\V-EiH\\o<Ch\ 

(5.186) 



(5.18c) 

||Hi - Hi;,||o < 

\\V-HiH\\o<Ch\ 

(5.18d) 


where E,;,, E^;,, Hih are the finite element solutions, k is the order of complete 
polynomials in the equal-order finite element interpolation. 


5.3 Time-Harmonic TE Waves 

For time-harmonic TE waves the first-order Maxwell’s equations sire 

dEy dEjc 


jwpi,Hz + 


_:^ = 0 inft. 


dx dy 


dH 

ju)e*E^ ^ = 0 in ft, 

ay 

dH 

j(jj£*Ey + = 0 in ft. 


dx 


dE^ . dE, 


+ 


V _ 


dx dy 


= 0 


in ft, 


(5.19a) 

(5.196) 

(5.19c) 


(5.19d) 
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in which e* = £r + jsi = £ + is the complex permittivity where the subscripts 

r and i indicate the real part and the imaginary part, respectively. For a complete 
description of TE wave problems appropriate boundary conditions should be included. 
One may consider, for example 

Hz = constant on F, (5.19e) 

+ Eyny = 0 on r, (5.19/) 

where n = (nx,ny) is the unit vector normal to the boundary F. The condition (5.19e) 
is an inhomogeneous version corresponding to (5.13g), and (5.19f) is a 2D version of 
(5.13h). We also remark that the boundary conditions (5.19e4) satisfy the boundary 
compatibility condition 

ju;e*(E,n* + Eyny) = ~ ^ (5.20) 

which is obtained by taking the operation n- to Elq. (5.19b) and (5.19c). 

In System (5.19) there are three unknowns and four equations, and thus the 
divergence-free equation (5.19d) seems redundant. By introducing a dummy variable 
i? into System (5.19) we have 


. „ ,dEy 

+ -gj- - 

= 0 

in n, 

(5.21a) 

. at? dHz 

+ ^ ay 

= 0 

in f2, 

(5.216) 

. W dH, 

j^cE, + ^+ 

= 0 

in Q, 

(5.21c) 

dEx dEy _ 
dx dy 

in 

fi. 

(5.21d) 


By taking the operation djdx to Eq. (5.21b) and the operation dJOy to Eq. (5.21c), 
and by adding the results together we obtain the Laplace equation for t?: 


dH dH 
dx^ dy^ 


in D. 


(5.22o) 


By taking the operation n- to the equations (5.21b) and (5.21c) and using the bound- 
ary compatibility condition (5.20) we obtain 


^ = 0 on F. (5.226) 

an 

From (5.22) we know that i? = constant^ that is, System (5.19) is completely equiv- 
alent to the augmented system (5.21) with four unknowns and four equations. Since 
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System (5.21) consists of two two-dimensional div-curl systems, and thus is elliptic. 
Therefore, System (5.19) is not ‘overdetermined’, but is indeed properly determined 
and dliptic. 

For numerical calctdation sepstrating the real and imaginary parts in (5.19a-d) 
leads to 


dEyr dE„ ^ 

dy 

(5.23a) 

Q CT 

u{erE^ + CiE„) ® 

(5.235) 

dH 

-w{erEyi -F CiEyr) + ^ 

(5.23c) 

dE„ , dEyr n . r> 

-H = 0 in Si, 

ox ay 

(5.23d) 

IT 1 dEyi dE,ei rt ' O 

(vuHrr + dy 

(5.23e) 

uj{trE„ - SiE^) - ^ = 0 

(5.23/) 

U:{£rEyr - SiEyi) + ^ = 0 171 fl. 

(5.23p) 

dExi dEyi . ^ 

dx dy 

(5.23h) 

course, in (5.23) the medium properties are different for different medium regions. 
We may write System (5.23) in the standaurd matrix form: 


in which 


9u 

Ai— -f- ■^2'^ AoU = f. 


/O 0 1 0 0 0\ 

0 0 0 0 0 0 

1 0 0 0 0 0 

0 1 0 0 0 0 

0 0 0 0 0 1 ’ 

0 0 0 0 0 0 

0 0 0 1 0 0 

^0 0 0 0 1 0 ; 


(5.24) 


(5.25a) 
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/ 0 -1 0 0 0 0 \ 

-1 0 0 0 0 0 

0 0 0 0 0 0 

0 010 00 

“ 0 0 0 0 -1 0 ’ 

0 0 0 -1 0 0 

0 0 0 0 0 0 

0 0 0 0 0 1 > 


0 

0 

0 

— W/i 

0 

0 \ 

0 

— we,' 

0 

0 

-we,. 

0 

0 

0 

— W€i 

0 

0 

—we. 

0 

0 

0 

0 

0 

0 

Ufl 

0 

0 

0 

0 

0 

0 

WCr 

0 

0 

— WCi 

0 

0 

0 


0 

0 

-WCi 

0 

0 

0 

0 

0 

0 / 


(5.256) 


(5.25c) 



At the interface Tint between two contiguous media (+) and (-) the following 
general conditions should be satisfied: 


nxE‘*'=nxE on Fin*, 


(5.26a) 


n X H+ = n X H- on Fi„*, (5.266) 

n • (e*+E+) = n • (e*“E”) on Fi„t, (5.26c) 

n . (/X+H+) = n • ifi-n-) on Tiru- (5-26d) 

For two-dimensional TE waves the interface conditions (5.26a) and (5.26c) become 

n,E+ - riyEt = n^E- - n^E' , (5.27a) 

n,E^ - UyE^ = n*E^ - nyE^, (5.276) 

n,(c+E+ - etEt) n,(e+E+ - efE^) = 


n^(e;E^ - €i E^) -|- n„(e, E^ - E^), 

n^Et + + nyiefE^ + e+Ej) = 


(5.27c) 
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^xr + ^xi) + ”v(^» ^yr + ^yi)t (5.27<f) 

The interface condition (5.26d) is automatically satisfied, and (5.26h) becomes 

H+ = H~. (5.27e) 

In the LSFEM the treatment of the interface conditions is not difficult. As in 
other node-based finite element methods, the nodes on the interface should be double- 
numbered. If a direct solver is employed for solving the discretized system, two ap- 
proaches are available: a simple way is to include the interface conditions into the 
least-squares functional; a better way is to use the conditions (5.27a-e) to modify the 
global stiffness matrix in the discretized system. If the conjugate gradient method 
is used, one just simply chooses the unknowns related to the medium (-f) (or (— )) 
as the true unknowns and keeps the conditions (5.27a-e) satisfied for each solution 
vector. 

Since the general formulation of the LSFEM has been given in Section 5.3, it is 
not necessary to write down the special one for the problem discussed in this section. 
One only needs to substitute the coefficients of (5.24) and the boundary conditions 
into a general-purpose LSFEM code. 

We consider two test problems that are taken from Paulsen and Lynch [49] in 
which the spurious solutions given by the curl-curl formulation as well as the correct 
solutions are illustrated. 

The first example is a cylinder (R=25) which is split into two regions having 
different complex permittivity. For the top re^on, = 3.0, e+ = —5.0 and = 1.0; 
for the bottom region, = 1.0, ej = 0.0 and fi~ = 1.0. This cylinder is excited 
by a uniform i?z|r = (1>0) (with u> = 0.05) imposed on the outer boundary. This 
problem is discretized by 932 bi-linear elements with 1016 nodes shown in Fig. 1(a). 
The contours of the computed real and imaginary magnetic field intensity are shown 
in Fig. 1(b) and (c), respectively. The vector plots of the reed and imaginary electric 
field intensity are illustrated in Fig. 1(d) and (e), respectively. 

In the second example, a smaller off-center cylinder (R=0.1) is embedded in a 
larger cylinder (R=0.25). The material properties for the outer region are e;!" = 0.0981, 
ef = —0.0196 and = 1.0; for the inner region e~ — 1-0, ej = 0.0 and = 1.0. 
A uniform ^*|r = (1,— 0.15) (with u = 44.7) is imposed on the outer boundary. 
Fig. 2(a) shows the mesh with 2027 bi-linear elements and 2165. The contours of 
the computed real and imaginary magnetic field intensity are shown in Fig. 2(b,c), 
respectively. The vector plots of the real and ima^nary electric field intensity are 
illustrated in Fig. l(d,e), respectively. 

As expected all computed results by the LSFEM are free of spurious modes. 

6 Conclusions 

(1) The system of the first-order full Maxwell’s equations seems “overspedfied” 
because it has more equations than unknowns. By taking into account of its div- 
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curl structure and introducing the dummy variables it proves to be properly deter- 
mined and elliptic in the space domain. The information provided by the divergence 
equations is not completely contained in the curl equations. Therefore, the diver- 
gence equations must be explicitly included in the first-order system to assure the 
uniqueness of the solution in steady-state cases and to guarantee the accuracy of the 
numerical solution for time-varying cases, 

(2) The least-squares method and the div-curl method are mathematically rigor- 
ous and useful tools for the derivation of correct second-order Maxwell’s equations 
and their boundary conditions. The curl-curl equations cannot stand alone, they must 
be supplemented by the divergence equations and additional natural boundary con- 
ditions to eliminate the spurious solutions. 

(3) The Helmholtz-type equations with appropriate natural boundary conditions, 
derived by the div-curl method or the least-squares method, can guarantee the im- 
plicit satisfaction of the divergence equations. For the solution of the Helmholtz-type 
equations the divergence conditions of the electric field and the magnetic field need to 
be enforced only on the electric wall and the magnetic symmetry wall, respectively. 

(4) The variational formulation corresponding to the Hdmholtz-type equations 
can be derived by using the least-squares look-alike method. This formulation theo- 
retically justifies that the penalty parameter in the Galerkin/penalty method should 
be taken as one. The advantage of this formulation is that the trial and test functions 
need only to satisfy the conditions related to the essential boundary conditions. 

(5) The node-based least-squares finite element method(LSFEM) can be used 
to solve both static and time-varying first-order Maxwdl’s equations directly and 
eflBciently with the divergence equations satisfied easily. The introduction of potentials 
and the gauging method, the edge element method, the staggered grid cuid upwinding, 
etc. all turn out to be unnecessary. 
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Fig. 1(a) The split cylinder and the mesh. 
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Fig. 1(b) Contours of constant magnetic field intensity Hr. 
Fig. 1(c) Contours of constant magnetic field intensity Hi. 
Fig. 1(d) Vectors of the computed electric field intensity E,. 
Fig. 1(e) Vectors of the computed electric field intensity E,. 
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Fig. 2(a) The off-center cylinder. 
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(d) (e) 


Fig. 2(b) Contours of constant magnetic field intensity Hr- 
Fig. 2(c) Contours of constant magnetic field intensity Hi- 
Fig. 2(d) Vectors of the computed electric field intensity E,.. 
Fig. 2(e) Vectors of the computed electric field intensity E*. 
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